Yellow eel length at age model summary

Load libraries

# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools","viridis","sf","patchwork","coda","boot","tidybayes","bayesplot", "nimble", "here","reshape2")

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  }

invisible(lapply(pkgs, library, character.only = T))

options(ggplot2.continuous.colour = "viridis")
theme_set(theme_minimal()) 

# Set path
home <- here::here()

Read data

## Read data
df.eel <- readRDS(file = "~/gitProjects/2024_R_DIASPARA/data/indhdb_3.RData") %>%
  st_drop_geometry()

Length at age - Schnute growth model

Define data

data.schnu <- df.eel %>%
  filter(fi_lfs_code == "Y") %>%
  rename(age = ageyear, 
         length = lengthmm) %>%
  drop_na(length, main_bas, age, ser_hty_code,ecoreg) %>%
  mutate(sex = as.integer(if_else(sex == "f", 0, 1)),
         emu = as.integer(factor(ser_emu_nameshort)),
         main_bas = as.integer(factor(main_bas)),
         er = as.integer(factor(ecoreg)),
         hab = as.integer(factor(ser_hty_code)),
         temp.sc = (tmp_dc_uyr - mean(tmp_dc_uyr))/sd(tmp_dc_uyr)) 

filter: removed 455,897 rows (52%), 416,562 rows remaining
rename: renamed 2 variables (length, age)
drop_na: removed 392,859 rows (94%), 23,703 rows remaining
mutate: converted 'sex' from character to integer (0 new NA)
        converted 'main_bas' from double to integer (0 new NA)
        new variable 'emu' (integer) with 38 unique values and 0% NA
        new variable 'er' (integer) with 10 unique values and 0% NA
        new variable 'hab' (integer) with 4 unique values and 0% NA
        new variable 'temp.sc' (double) with 66 unique values and 0% NA

Plot length at age data

# density of length by sex
data.schnu %>%
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot() +
  geom_density(aes(x = age, fill = sex), alpha = 0.3) +
  theme_minimal() +
  
data.schnu %>% 
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot() +
  geom_density(aes(x = length, fill = sex), alpha = 0.3) +
  theme_minimal()

mutate: converted 'sex' from integer to character (0 new NA)
mutate: converted 'sex' from integer to character (0 new NA)

# length at age by sex
data.schnu %>%
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot() +
  geom_point(aes(age, length, color = factor(sex)), size = 0.1, alpha = 0.8) +
  facet_wrap(~ecoreg) +
  theme_minimal()

mutate: converted 'sex' from integer to character (0 new NA)

# yellow growth by habitat and region, lots of combinations missing 
data.schnu %>%
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot() +
  geom_histogram(aes(length), size = 0.1, alpha = 0.8) +
  facet_grid(ser_hty_code ~ ecoreg, scales = "free_y") +
  scale_x_continuous(breaks = seq(0, 1000, 500)) +
  theme_minimal() 

mutate: converted 'sex' from integer to character (0 new NA)

Warning in geom_histogram(aes(length), size = 0.1, alpha = 0.8): Ignoring
unknown parameters: `size`

`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# individuals (obs) per basin
data.schnu %>%
  summarise(nobs.bas = n(), .by = c(main_bas, ecoreg)) %>%
  ggplot() +
  geom_boxplot(aes(x = ecoreg, y = nobs.bas), show.legend = FALSE) +
  geom_point(aes(ecoreg, nobs.bas), alpha = 0.5, show.legend = FALSE) +
  theme_minimal() +
  labs(caption = "observations per basin by region") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = rel(.75)))

summarise: now 122 rows and 3 columns, ungrouped

# possible temperature effects on length at age??
data.schnu %>%
  ggplot(aes(tmp_dc_uyr/10, length)) +
  geom_point(show.legend = FALSE, alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE) +
  facet_wrap(~factor(round(age))) +
  theme_minimal() +
  labs(x = "Temp. [C]")

`geom_smooth()` using formula = 'y ~ x'

Warning in qt((1 - level)/2, df): NaNs produced

Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf

Source Schnute model / Load samples

# Load samples 
schnu.samples <- readRDS(paste0(home,"/models_eel/samples/schnu12_samples_2026-01-23.RData"))$samples

# # Or source model code and mcmc (load libs and data)
# t <- Sys.time()
# source(file = paste0(home,"/models_eel/model_eel_lengthatage.R"))
# Sys.time() - t # approx 15 hours with 10000 hmc smaples

Traces and ‘potential scale reduction factor’

The variables in the Schnute model are: “ps”,“l.sig”,“Ustar”,“U”,“hK”,“mu.par”,“sig.par”,“btL”,“btK”,“sex”,“bsL”,“bsK”,“par”,“l.mu”),

## par ( a sample of the par traces)
node.sub <- grep("^par\\[", colnames(schnu.samples[[1]]), value = TRUE)
schnu.samples[, node.sub[sample(1:length(node.sub), 36)], drop = FALSE] %>%
mcmc_trace()

pl <- schnu.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
# which psrf are > 1.01 
pl$psrf[which(pl$psrf[, 1] > 1.01),]

        Point est. Upper C.I.
par[4]    1.048458   1.137373
par[8]    1.011478   1.011845
par[13]   1.045924   1.149102
par[16]   1.016553   1.040680
par[18]   1.021420   1.047024
par[19]   1.012539   1.058384
par[20]   1.022163   1.056304
par[23]   1.023025   1.089679
par[25]   1.011919   1.012801
par[26]   1.017843   1.030151
par[27]   1.038023   1.072037
par[28]   1.020924   1.022187
par[32]   1.032481   1.032553
par[33]   1.018023   1.047235
par[36]   1.021039   1.022934
par[39]   1.057881   1.239196
par[40]   1.012711   1.020552

# check traces of those psrf that are > 1.01
schnu.samples[,node.sub[which(pl$psrf[, 1] > 1.01)]] %>%
#schnu.samples[, node.sub[sample(1:length(node.sub), 16)], drop = FALSE] %>%
mcmc_trace()

# EES of par 
pl <- schnu.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
pl[which(pl < 1000)]

  par[1]  par[11]  par[12]  par[16]  par[17]  par[19]  par[20]  par[31] 
843.7755 570.4663 761.8103 972.5341 942.5795 552.4972 468.0654 689.9927 
 par[32]  par[39]  par[40] 
422.0575 392.1275 352.9016 

## bpar ( a sample of the par traces)
node.sub <- grep("^bpar\\[", colnames(schnu.samples[[1]]), value = TRUE)
schnu.samples[, node.sub[sample(1:length(node.sub), 36)], drop = FALSE] %>%
mcmc_trace()

pl <- schnu.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
# which psrf are > 1.01 
pl$psrf[which(pl$psrf[, 1] > 1.1),] # bpar[8, 3] above 1.1

Point est. Upper C.I. 
  1.147959   1.350098 

# check traces of those psrf that are > 1.01
schnu.samples[,node.sub[which(pl$psrf[, 1] > 1.1)]] %>%
#schnu.samples[, node.sub[sample(1:length(node.sub), 16)], drop = FALSE] %>%
mcmc_trace()

# EES of par 
pl <- schnu.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
pl[which(pl < 1000)]

  bpar[3, 1]   bpar[5, 1]   bpar[6, 1]   bpar[7, 1]   bpar[8, 1]  bpar[14, 1] 
    948.1432     770.9566     950.5805     857.6619     753.4540     946.3053 
 bpar[15, 1]  bpar[17, 1]  bpar[22, 1]  bpar[84, 1]  bpar[87, 1] bpar[122, 1] 
    676.4296     895.4444     838.2759     727.3626     857.0061     943.3211 
  bpar[1, 2]   bpar[2, 2]   bpar[3, 2]   bpar[4, 2]   bpar[5, 2]   bpar[6, 2] 
    601.3093     741.2412     700.8293     750.6786     429.3968     601.3795 
  bpar[7, 2]   bpar[8, 2]   bpar[9, 2]  bpar[10, 2]  bpar[11, 2]  bpar[12, 2] 
    602.9181     872.6277     863.0809     919.1002     870.6416     845.8967 
 bpar[13, 2]  bpar[14, 2]  bpar[15, 2]  bpar[16, 2]  bpar[17, 2]  bpar[19, 2] 
    690.0612     849.2242     761.9672     977.0764     654.7935     995.5587 
 bpar[22, 2]  bpar[32, 2]  bpar[35, 2]  bpar[36, 2]  bpar[39, 2]  bpar[40, 2] 
    850.5284     916.3299     917.0442     720.3041     817.2068     633.9661 
 bpar[42, 2]  bpar[43, 2]  bpar[44, 2]  bpar[45, 2]  bpar[46, 2]  bpar[47, 2] 
    889.4795     938.3641     893.7295     957.1515     955.1212     648.8139 
 bpar[48, 2]  bpar[49, 2]  bpar[51, 2]  bpar[52, 2]  bpar[54, 2]  bpar[55, 2] 
    786.0697     824.0975     745.0648     848.4683     731.7797     976.7650 
 bpar[56, 2]  bpar[57, 2]  bpar[58, 2]  bpar[64, 2]  bpar[66, 2]  bpar[78, 2] 
    779.4927     935.0357     926.5370     860.6092     991.1361     848.4313 
 bpar[79, 2]  bpar[80, 2]  bpar[84, 2]  bpar[86, 2]  bpar[87, 2]  bpar[88, 2] 
    901.5305     919.3909     817.8378     980.6844     659.0626     872.0237 
 bpar[91, 2]  bpar[92, 2]  bpar[94, 2]  bpar[95, 2]  bpar[96, 2]  bpar[97, 2] 
    486.1416     669.4739     803.7739     555.7704     649.1016     755.1929 
 bpar[98, 2] bpar[100, 2] bpar[102, 2] bpar[103, 2] bpar[105, 2] bpar[109, 2] 
    428.6735     851.0072     644.4590     825.9560     978.6525     873.0965 
bpar[110, 2] bpar[114, 2] bpar[120, 2]  bpar[39, 3]  bpar[64, 3]  bpar[91, 3] 
    738.6507     941.6395     763.0571     932.2886     963.1855     893.0814 
 bpar[98, 3]   bpar[1, 4]   bpar[3, 4]   bpar[4, 4]   bpar[5, 4]   bpar[6, 4] 
    883.9922     820.7666     999.8462     858.0152     698.5759     800.9391 
  bpar[7, 4]   bpar[8, 4]   bpar[9, 4]  bpar[10, 4]  bpar[11, 4]  bpar[12, 4] 
    812.7339     818.6262     883.4761     991.4536     856.7663     984.8809 
 bpar[13, 4]  bpar[15, 4]  bpar[16, 4]  bpar[17, 4]  bpar[19, 4]  bpar[32, 4] 
    879.7664     810.7132     925.2374     867.9584     865.6049     885.7319 
 bpar[36, 4]  bpar[37, 4]  bpar[38, 4]  bpar[39, 4]  bpar[40, 4]  bpar[41, 4] 
    636.3827     799.0061     767.6049     460.4679     643.1860     814.0988 
 bpar[42, 4]  bpar[43, 4]  bpar[44, 4]  bpar[45, 4]  bpar[46, 4]  bpar[47, 4] 
    603.5091     856.2787     670.7697     641.8752     652.8818     442.2811 
 bpar[48, 4]  bpar[49, 4]  bpar[50, 4]  bpar[51, 4]  bpar[52, 4]  bpar[53, 4] 
    500.0390     894.6545     675.7451     607.1579     659.1110     652.1654 
 bpar[54, 4]  bpar[55, 4]  bpar[56, 4]  bpar[57, 4]  bpar[58, 4]  bpar[59, 4] 
    733.4089     745.4376     639.3637     543.8596     463.1725     567.4424 
 bpar[60, 4]  bpar[61, 4]  bpar[62, 4]  bpar[63, 4]  bpar[64, 4]  bpar[65, 4] 
    844.9550     660.7177     952.0253     929.4329     807.6663     816.9976 
 bpar[66, 4]  bpar[67, 4]  bpar[68, 4]  bpar[69, 4]  bpar[73, 4]  bpar[74, 4] 
    893.4789     643.0557     817.4885     650.9028     580.8305     798.2366 
 bpar[75, 4]  bpar[76, 4]  bpar[77, 4]  bpar[78, 4]  bpar[79, 4]  bpar[80, 4] 
    887.6696     733.9682     913.7768     810.0024     844.6848     595.5053 
 bpar[81, 4]  bpar[86, 4]  bpar[87, 4]  bpar[88, 4]  bpar[89, 4]  bpar[90, 4] 
    642.5883     687.0014     480.4702     386.3373     529.4692     529.8801 
 bpar[91, 4]  bpar[92, 4]  bpar[93, 4]  bpar[94, 4]  bpar[95, 4]  bpar[96, 4] 
    424.9476     487.7709     564.4911     530.1805     521.7596     517.4336 
 bpar[97, 4]  bpar[98, 4]  bpar[99, 4] bpar[100, 4] bpar[101, 4] bpar[102, 4] 
    501.0609     424.1930     475.9285     517.8181     601.6196     566.6986 
bpar[103, 4] bpar[104, 4] bpar[105, 4] bpar[106, 4] bpar[107, 4] bpar[108, 4] 
    589.9626     535.8460     553.8243     555.7366     554.1357     655.4609 
bpar[109, 4] bpar[110, 4] bpar[111, 4] bpar[112, 4] bpar[120, 4] bpar[121, 4] 
    483.3110     488.9217     555.2188     523.8316     975.0047     619.3883 

## mu.par
node.sub <- grep("^mu.par", colnames(schnu.samples[[1]]), value = TRUE)
schnu.samples[, node.sub[], drop = FALSE] %>%
mcmc_trace()

pl <- schnu.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
# no mu.par psrf > 1.01
pl$psrf[which(pl$psrf[, 1] > 1.01),]

           Point est. Upper C.I.
mu.par[1]    1.011881   1.058185
mu.par[2]    1.035215   1.069764
mu.par[4]    1.018474   1.062684
mu.par[18]   1.011681   1.012237
mu.par[26]   1.025403   1.081335
mu.par[35]   1.026952   1.072670
mu.par[39]   1.018520   1.086508

# all ess are good
schnu.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()

 mu.par[1]  mu.par[2]  mu.par[3]  mu.par[4]  mu.par[5]  mu.par[6]  mu.par[7] 
  909.4438  1242.6062  1808.1295  1613.8457  2441.1963  1960.5424  2078.5496 
 mu.par[8]  mu.par[9] mu.par[10] mu.par[11] mu.par[12] mu.par[13] mu.par[14] 
 2442.6133  1667.3980  2343.3669  1640.0217  2195.5348  1992.8654  2427.8127 
mu.par[15] mu.par[16] mu.par[17] mu.par[18] mu.par[19] mu.par[20] mu.par[21] 
 2749.5633  2762.5897  1612.5406  2514.5778  2871.8270  2695.4313  1876.0887 
mu.par[22] mu.par[23] mu.par[24] mu.par[25] mu.par[26] mu.par[27] mu.par[28] 
 2550.5469  1833.1477  2098.3114  2492.4622  2371.2556  2120.4373  2564.0103 
mu.par[29] mu.par[30] mu.par[31] mu.par[32] mu.par[33] mu.par[34] mu.par[35] 
 1575.4385  2395.8969   951.8449  2373.0690  2309.8653  2071.4054  1546.5127 
mu.par[36] mu.par[37] mu.par[38] mu.par[39] mu.par[40] 
 2146.4663  1563.1707  2218.7569   532.2122  1302.4909 

## sig.par
node.sub <- grep("^sig.par", colnames(schnu.samples[[1]]), value = TRUE)
schnu.samples[, node.sub[], drop = FALSE] %>%
mcmc_trace()

schnu.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()

Potential scale reduction factors:

            Point est. Upper C.I.
sig.par[1]        1.01       1.05
sig.par[2]        1.04       1.12
sig.par[3]        1.01       1.02
sig.par[4]        1.00       1.00
sig.par[5]        1.01       1.04
sig.par[6]        1.00       1.03
sig.par[7]        1.00       1.01
sig.par[8]        1.00       1.02
sig.par[9]        1.01       1.04
sig.par[10]       1.02       1.08
sig.par[11]       1.01       1.01
sig.par[12]       1.00       1.00
sig.par[13]       1.01       1.05
sig.par[14]       1.01       1.01
sig.par[15]       1.00       1.00
sig.par[16]       1.00       1.01
sig.par[17]       1.00       1.01
sig.par[18]       1.00       1.02
sig.par[19]       1.00       1.00
sig.par[20]       1.02       1.05
sig.par[21]       1.00       1.01
sig.par[22]       1.02       1.09
sig.par[23]       1.00       1.00
sig.par[24]       1.00       1.00
sig.par[25]       1.01       1.02
sig.par[26]       1.03       1.11
sig.par[27]       1.01       1.03
sig.par[28]       1.04       1.09
sig.par[29]       1.00       1.00
sig.par[30]       1.00       1.01
sig.par[31]       1.00       1.01
sig.par[32]       1.02       1.09
sig.par[33]       1.02       1.07
sig.par[34]       1.03       1.11
sig.par[35]       1.00       1.00
sig.par[36]       1.03       1.11
sig.par[37]       1.02       1.11
sig.par[38]       1.00       1.01
sig.par[39]       1.00       1.02
sig.par[40]       1.04       1.19

Multivariate psrf

1.15

schnu.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()

 sig.par[1]  sig.par[2]  sig.par[3]  sig.par[4]  sig.par[5]  sig.par[6] 
   365.6773    531.1351   1826.2333   1874.6235   2553.7545   2307.2927 
 sig.par[7]  sig.par[8]  sig.par[9] sig.par[10] sig.par[11] sig.par[12] 
  1361.0411   2016.8317    848.0077   1297.8696    337.6132   1250.2397 
sig.par[13] sig.par[14] sig.par[15] sig.par[16] sig.par[17] sig.par[18] 
  1271.4862   1891.8112   1719.6556   1306.0756   1270.3260   1873.6115 
sig.par[19] sig.par[20] sig.par[21] sig.par[22] sig.par[23] sig.par[24] 
   880.1053    900.9995   1990.6687   1351.5744   1427.8207   1747.3307 
sig.par[25] sig.par[26] sig.par[27] sig.par[28] sig.par[29] sig.par[30] 
  1645.1808   1434.1062   1456.1288   2208.7470   2230.2276    974.8125 
sig.par[31] sig.par[32] sig.par[33] sig.par[34] sig.par[35] sig.par[36] 
   686.3089    422.0669   1366.0397   1993.7411   1540.8427    873.7784 
sig.par[37] sig.par[38] sig.par[39] sig.par[40] 
  1356.1762   2343.7182    321.3134    518.7500 

# R - cholesky of correlation matrix
node.sub <- grep("^R\\[", colnames(schnu.samples[[1]]), value = TRUE)
schnu.samples[, node.sub[sample(1:length(node.sub), 48)], drop = FALSE] %>%
mcmc_trace()

pl <- schnu.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
pl[which(pl < 1000 & pl > 0)]

named numeric(0)

# l.sig bsL bsK btK btL
node.sub <- grep("^(l.sig|bsL|bsK|btK|btL)", colnames(schnu.samples[[1]]), value = TRUE)
schnu.samples[, node.sub[], drop = FALSE] %>%
mcmc_trace()

schnu.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()

Potential scale reduction factors:

      Point est. Upper C.I.
bsK         1.00       1.02
bsL         1.01       1.01
btK         1.00       1.00
btL         1.00       1.00
l.sig       1.00       1.01

Multivariate psrf

1.01

schnu.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()

      bsK       bsL       btK       btL     l.sig 
 403.8104 3000.0000  683.1211 2778.7660 2759.0845 

# ps
node.sub <- grep("^ps", colnames(schnu.samples[[1]]), value = TRUE)
schnu.samples[, node.sub[], drop = FALSE] %>%
mcmc_trace()

schnu.samples[, node.sub[], drop = FALSE] %>%
gelman.diag(multivariate = FALSE)

Potential scale reduction factors:

      Point est. Upper C.I.
ps[1]       1.02       1.07
ps[2]       1.02       1.07

schnu.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()

   ps[1]    ps[2] 
1704.048 1704.048 

node.sub <- grep("^hK", colnames(schnu.samples$chain1), value = TRUE)
schnu.samples[, node.sub[], drop = FALSE] %>%
mcmc_trace()

pl <- schnu.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
pl$psrf[which(pl$psrf[, 1] > 1.01),]

      Point est. Upper C.I.
hK[1]   1.039077   1.170452
hK[2]   1.014245   1.068827
hK[3]   1.074550   1.296397
hK[4]   1.011147   1.055427

Schnute pars and covariates

Schnute pars

dfp <- schnu.samples %>%
  gather_draws(bpar[main_bas,pars], sep = ",") %>%
  mutate(par = case_when(pars == 1 ~ "L1",
                         pars == 2 ~ "P",
                         pars == 3 ~ "L2",
                         pars == 4 ~ "K"),
         val = case_when(!pars == 2 ~ exp(.value),
                         .default = .value)) %>%
  left_join(data.schnu %>% distinct(ecoreg, main_bas))

mutate (grouped): new variable 'par' (character) with 4 unique values and 0% NA
                  new variable 'val' (double) with 1,462,048 unique values and 0% NA
distinct: removed 23,581 rows (99%), 122 rows remaining
Joining with `by = join_by(main_bas)`
left_join: added one column (ecoreg)
           > rows only in x                                  0
           > rows only in data.schnu %>% distinct.. (        0)
           > matched rows                            1,464,000
           >                                        ===========
           > rows total                              1,464,000

dfp %>%
  filter(par == "L1") %>%
  ggplot() +
  geom_density(aes(x = val, group = main_bas), color = "grey") + 
  geom_density(aes(x = val), color = "black", show.legend = FALSE) + 
  facet_wrap(~ecoreg, scales = "free", ncol = 3) +
  labs(title = "L1 by main basin in region", caption = "black line is region specific posterior")

filter (grouped): removed 1,098,000 rows (75%), 366,000 rows remaining (removed
366 groups, 122 groups remaining)

ggsave(paste0(home, "/report/elaa_L1.png"), width = 17, height = 14, units = "cm")

dfp %>%
  filter(par == "L2") %>%
  ggplot() +
  geom_density(aes(x = val, group = main_bas), color = "grey") + 
  geom_density(aes(x = val), color = "black", show.legend = FALSE) + 
  facet_wrap(~ecoreg, scales = "free", ncol = 3) +
  labs(title = "L2 by main basin in region", caption = "black line is region specific posterior")

filter (grouped): removed 1,098,000 rows (75%), 366,000 rows remaining (removed
366 groups, 122 groups remaining)

ggsave(paste0(home, "/report/elaa_L2.png"), width = 17, height = 14, units = "cm")

dfp %>%
  filter(par == "K") %>%
  ggplot() +
  geom_density(aes(x = val, group = main_bas), color = "grey") + 
  geom_density(aes(x = val), color = "black", show.legend = FALSE) + 
  facet_wrap(~ecoreg, scales = "free", ncol = 3) +
  labs(title = "K by main basin in region", caption = "black line is region specific posterior")

filter (grouped): removed 1,098,000 rows (75%), 366,000 rows remaining (removed
366 groups, 122 groups remaining)

ggsave(paste0(home, "/report/elaa_K.png"), width = 17, height = 14, units = "cm")

dfp %>%
  filter(par == "P") %>%
  ggplot() +
  geom_density(aes(x = val, group = main_bas), color = "grey") + 
  geom_density(aes(x = val), color = "black", show.legend = FALSE) + 
  facet_wrap(~ecoreg, scales = "free", ncol = 3) +
  labs(title = "P by main basin in region", caption = "black line is region specific posterior")

filter (grouped): removed 1,098,000 rows (75%), 366,000 rows remaining (removed
366 groups, 122 groups remaining)

ggsave(paste0(home, "/report/elaa_P.png"), width = 17, height = 14, units = "cm")

Covariates

schnu.samples %>%
  gather_draws(hK[hab], sep = ",") %>%
  ungroup() %>%
  left_join(data.schnu %>% distinct(hab,ser_hty_code)) %>%
  #filter(!is.na(ecoreg)) %>%
  ggplot() +
  geom_density(aes(x = .value, color = factor(ser_hty_code))) +
  #facet_wrap(~ecoreg, scales = "free", ncol = 3) +
  labs(title = "habitat")

ungroup: no grouping variables remain
distinct: removed 23,699 rows (>99%), 4 rows remaining
Joining with `by = join_by(hab)`
left_join: added one column (ser_hty_code)
           > rows only in x                               0
           > rows only in data.schnu %>% distinct.. (     0)
           > matched rows                            12,000
           >                                        ========
           > rows total                              12,000

schnu.samples %>%
  gather_draws(btL,btK,bsL,bsK) %>%
  ggplot() +
  geom_density(aes(x = .value, color = .variable)) +
  facet_wrap(~.variable, scales = "free")

ggsave(paste0(home, "/report/elaa_bLK.png"), width = 17, height = 14, units = "cm")

Predicted / observed

dfp <- schnu.samples %>%
  spread_draws(l.mu[i]) %>%
  median_qi() %>%
  rename(pred.l = l.mu,
         ind.id = i) %>%
  ungroup() %>%
  left_join(data.schnu %>% 
              select(length,ecoreg) %>% 
              rename(obs.l = length) %>%
              mutate(ind.id = row_number())) 

rename: renamed 2 variables (ind.id, pred.l)
ungroup: no grouping variables remain
select: dropped 348 variables (mei_fi_id, fi_year, fi_date, ser_x, ser_y, …)
rename: renamed one variable (obs.l)
mutate: new variable 'ind.id' (integer) with 23,703 unique values and 0% NA
Joining with `by = join_by(ind.id)`
left_join: added 2 columns (obs.l, ecoreg)
           > rows only in x                               0
           > rows only in data.schnu %>% select(l.. (     0)
           > matched rows                            23,703
           >                                        ========
           > rows total                              23,703

dfp %>%
  slice_sample(n = 5000) %>%
  ggplot() +
  geom_abline(slope =1, color = "darkgrey") +
  geom_point(aes(x = obs.l, y = pred.l), shape = 4, alpha = 0.7) +
  facet_wrap(~ecoreg) +
  theme_minimal() 

slice_sample: removed 18,703 rows (79%), 5,000 rows remaining

ggsave(paste0(home, "/report/elaa_po.png"), width = 17, height = 14, units = "cm")

Predicted Schnute curves

dfp <- schnu.samples[[1]] %>% # use one chain to reduce object size
  spread_draws(l.mu[i], sex[i]) %>%
  ungroup() %>%
  left_join(data.schnu %>% select(main_bas, ser_hty_code, ecoreg, age, length) %>% mutate(id = row_number()), join_by(i == id))

ungroup: no grouping variables remain
select: dropped 345 variables (mei_fi_id, fi_year, fi_date, ser_x, ser_y, …)
mutate: new variable 'id' (integer) with 23,703 unique values and 0% NA
left_join: added 5 columns (main_bas, ser_hty_code, ecoreg, age, length)
           > rows only in x                                   0
           > rows only in data.schnu %>% select(m.. (         0)
           > matched rows                            35,554,500
           >                                        ============
           > rows total                              35,554,500

dfp.summb <- dfp %>%
  summarise(medl = median(l.mu), .by = c(sex, age, ecoreg, main_bas, ser_hty_code))

summarise: now 1,618 rows and 6 columns, ungrouped

dfp %>%
  filter(sex == 0) %>%
  ungroup() %>%
  ggplot() +
  geom_boxplot(aes(x = age, y = l.mu, group = age), 
               outlier.color = "lightgrey",
               outlier.fill = "lightgrey",
               outlier.size = 0.5
               ) +
  geom_line(data = dfp.summb %>% filter(sex == 0), 
            aes(x = age, y = medl, group = factor(main_bas), color = "red"), alpha = 0.5, show.legend = FALSE) +
  facet_wrap(ser_hty_code~ecoreg) +
  theme_minimal() +
  labs(title = "females")

filter: removed 3,188,645 rows (9%), 32,365,855 rows remaining
ungroup: no grouping variables remain
filter: removed 418 rows (26%), 1,200 rows remaining

ggsave(paste0(home, "/report/elaa_fc.png"), width = 17, height = 14, units = "cm")

dfp %>%
  filter(sex == 1) %>%
  ungroup() %>%
  ggplot() +
  geom_boxplot(aes(x = age, y = l.mu, group = age),
               outlier.color = "lightgrey",
               outlier.fill = "lightgrey",
               outlier.size = 0.5
               ) +
  geom_line(data = dfp.summb %>% filter(sex == 1), 
            aes(x = age, y = medl, group = factor(main_bas), color = "red"), alpha = 0.5, show.legend = FALSE) +
  facet_wrap(ser_hty_code~ecoreg) +
  theme_minimal() +
  labs(title = "males")

filter: removed 32,365,855 rows (91%), 3,188,645 rows remaining
ungroup: no grouping variables remain
filter: removed 1,200 rows (74%), 418 rows remaining
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

ggsave(paste0(home, "/report/elaa_mc.png"), width = 17, height = 14, units = "cm")

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?